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We consider the fracture of a free-standing two-dimensional (2D) elastic-brittle network to be 
used as protective coating subject to constant tensile stress applied on its rim. Using a Molecular 
Dynamics simulation with Langevin thermostat, we investigate the scission and recombination of 
bonds, and the formation of cracks in the 2D graphene-like hexagonal sheet for different pulling 
force / and temperature T. We find that bond rupture occurs almost always at the sheet periphery 
and the First Mean Breakage Time (r) of bonds decays with membrane size as (t) oc where 
P « 0.50 ± 0.03 and A'^ denotes the number of atoms in the membrane. The probability distribution 

■ of bond scission times t is given by a Poisson function W{t) oc t^^'^ exp(— t/(r)). The mean failure 
time (rr) that takes to rip-off the sheet declines with growing size A'' as a power law (r^) oc N^'^'^^K 

■ We also find (r^) oc exp{AUo / ksT) where the nucleation barrier for crack formation AUo oc /^^, in 
agreement with Griffith's theory, (rr) displays an Arrhenian dependence of (r^) on temperature T. 
Our results indicate a rapid increase in crack spreading velocity with growing external tension /. 
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I. INTRODUCTION 

c/2 ■ Fracture in engineering materials is a long-standing topic of research due to problems that arise with technological 
applications and the ensuing economical implications. Thus, for decades a lot of attention has been focused on 
understanding the macroscopic and microscopic factors which trigger failure. Recently, the interest and the need for 
better understanding of the interplay between elastic and fracture properties of brittle materials has been revived due 
to the rapidly developing design of advanced structural materials. 

Promising aspects for applications include reversible polymer networks 0,07 and also graphene, that shows unusual 
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Q I thermomechanical properties [1, 01 • Among other things, graphene, which is a honey-comb lattice packed with C 
atoms can be used as anti-corrosion gas barrier protective coating Q, in chemical and bio-sensors @, or as efficient 
membrane for gas separation 0. In aU possible Applications the terTperature and stress-dependent fracture strength 
, of this 2D-network is of crucial importance. Graphene has been investigated recently by Barnard and Snook Q using 
^ ■ ab initio quantum mechanical techniques whereby it was noted that that the problems "has been overlooked by most 
^\ , computational and theoretical studies" . 

. An important example of biological microstructure is spectrin, the red blood cell membrane skeleton, which reinforces 
] the cytoplasmic face of the membrane. In erytrhrocytes, the membrane skeleton enables it to undergo large extensional 
deformations while maintaining the structural integrity of the membrane. A number of studies, based on continuum- 
0, percolation- [lol - [T2j . or molecular level [iq considerations of the mechanical breakdown of this network, 
modeled as a triangular lattice of spectrin tetramers, have been reported so far. Many of these studies can be viewed 
in a broader context as part of the problem of thermal decomposition of gels , epoxy resins [3 and other 3D 
networks both experimentally [l6l - [Tg |. and by means of simulations [l^ in the case of Poly-dimethylsiloxane (PDMS). 

The afore-mentioned examples illustrate well the need for deeper understanding of the processes of failure in brittle 
materials. Besides analytical and laboratory investigations, computer simulations (20l - [2^ have provided meanwhile a 
lot of insight in aspects that are difficult for direct observations or theoretical treatment - for a review of previous works 



X 



'—I 

d ' see Alava et al. [23[. Most of these studies focus on the propagation of (pre-existing) cracks, relating observations to 
the well known Griffith's model [24| of crack formation. A number of important aspects of material failure have found 
thereby little attention. Thus only a few simulations examine the rate of crack nucleation which involves long time 
scales necessary for thermal activation - see, however, [25l428| . Effects of system size on the characteristic time for bond 
rupture have not been examined except in a recent MD study by Dias et al. [29| . Also recombination of broken bonds 
has not been considered. These and other insufficiently explored properties related to fracture have motivated our 
present investigation of a free-standing 2D honeycomb brittle membrane by means of Molecular Dynamics simulation. 
In view of the possible applications as anti-corrosion and gas barrier coating, we consider a radially-spanned sheet 
of regular hexagonal flake shape so as to minimize effects of corners and unequal edge lengths that are typical for 
ribbon-like sheets. Tensile constant force is applied on the rim of the flake, perpendicular to each edge. By varying 
system size, tensile force and temperature, we collect a number of results which characterize the initiation and the 
course of fragmentation in stretched 2D honeycomb networks. 

The paper is organized as follows: after a brief introduction, we sketch our model in Sec. |TT] where we consider 
interactions between atoms in the brittle honeycomb membrane, define the threshold for bond scission, and also 
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introduce some basic quantities that are measured in the course of the simulation. In Sec. IIIII we present our 
simulation results, presenting briefly the results on recombination of broken bonds - IIII Al the distribution of bond 
scission rates over the membrane surface, the dependence of the Mean First Breakage Time (MFBT) before a bond 
scission takes place and of the mean failure time until the 2D sheet breaks apart on applied tensile force, and examine 
how these times depend on membrane size and temperature - IIIIBI The formation of cracks at different cases of 
applied stress as well as their propagation in a 2D honeycomb brittle sheet are briefly considered in subsection IIII CI 
We end this report by a brief summary of results in Section IIVI 

II. MODEL AND SIMULATION PROCEDURE 
A. The model 

We study a coarse-grained model of honeycomb membrane embedded in three-dimensional (3D) space. The mem- 
brane consists of N spherical particles (beads, monomers) of diameter a connected in a honey-comb lattice structure 
whereby each monomer is bonded with three nearest-neighbors except for the monomers on the membrane edges 
which have only two bonds (see Fig. [T] [left panel]). The total number of monomers N in such a membrane is = 
where by L we denote the number of monomers (or hexagon cells) on the edge of the membrane (i.e., L characterizes 
the linear size of the membrane). There are altogether Nbonds = {-^N — 6L)/2 bonds in the membrane. In our studies 
we consider symmetric hexagonal membranes (i.e., flakes) so as to minimize possible effects due to the asymmetric of 
edges or vortices at the membrane periphery. 




FIG. 1: [left panel] A membrane with honeycomb structure that contains a total of A'^ = 54 beads and has linear size L — 3 {L 
is the number of hexagonal cells on the edge of the membrane), [right panel] A snapshot of a typical conformation of an intact 
membrane with L = 30 containing 5400 monomers after equilibration with no external force applied. Typical wrinkles are seen 
to form on the surface. 

For the analysis of our results we find it appropriate to divide the two-dimensional membrane network so that all 
bonds fall into different subgroups presented by concentric "circles" with consecutive numbers (sec Fig. [T] [right panel]) 
proportional to their radial distance from the membrane center. To odd circle numbers thus belong bonds that are 
nearly tangential to the corresponding circle. Even circles contain no encompass radially oriented bonds (shown to 
cross the circle in Fig. [T]). The total number of circles C in a membrane of linear size L is found to be C = (2L — 1). 
We use this scheme of labeling the groups of bonds that compose the membrane in order to represent our simulation 
results in appropriate way which relates them to their relative proximity to membrane's periphery. 

B. Potentials 

The nearest-neighbors in the membrane are connected to each other by breakable anharmonic bonds described by 
a Morse potential, 

Uuir) = eA/{l - exp[-a(r - n„i„)]}^ (1) 

where r is the distance between the monomers. Here a = 1 is a constant that determines the width of the potential 
well (i.e., bond elasticity) and r^i^ = 1 is the equilibrium bond length. The dissociation energy of a given bond. 
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Em = 1, is measured in units of ksT where ks denotes the Boltzmann constant and T is the temperature. The 
minimum of this potential occurs at r = r^i^, C^Morsol^'min) = 0. The maximal restoring force of the Morse potential, 
/max = —dUM/dr = Q;e7v//2. is reached at the inflection point, r = ?'niin + Q;~^ lii(2) ~ 2.69. This force /max determines 
the maximal tensile strength of the membranes bonds. Since Um{0) ~ 2.95, the Morse potential, Eq. ([T]), is only 
weakly repulsive and beads could partially penetrate one another at r < Tmin- Therefore, in order to allow properly 
for the excluded volume interactions between bonded monomers, we take the bond potential as a sum of Uuir) and 
the so called Weeks-Chandler- Anderson (WCA) potential, C/wca(?'), (i-C-, the shifted and truncated repulsive branch 
of the Lennard- Jones potential). 
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with parameter e = 1 and monomer diameter a = w 0.89 so that the minimum of the WCA potential to 

coincides with the minimum of the Morse potential. Thus, the length scale is set by the parameter rmin = '2}l^a — 1. 
The nonbonded interactions between monomers are taken into account by means of the WCA potential, Eq. Thus, 
the nonbonded interactions in our model correspond to good solvent conditions whereas the bonded interactions make 
the bonds breakable when subject to stretching. External stretching force / is applied to monomers at the membrane 
rim in direction perpendicular to the respective edge - Fig. [2^. 

Before we turn to the problem of membrane failure under constant tensile force, we show here some typical elastic 
properties of the intact honeycomb network sheet that is used in our computer experiments - Fig. [2] In Fig. 
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FIG. 2: (a) A protective honeycomb network is spanned at the orifice of a prism whose size may vary due to thermal expansion. 
Tensile forces acting on the membrane periphery are indicated by red arrows, (b) Mean strain of a honeycomb membrane of 
size L = 10 as a function of external tensile stress / at T = 0.01 and 7 — 0.25. Depending on the way in which external 
force is applied: (i) radial stretching - c.f. Fig. [Of, or (ii) uniaxial stretching - Fig. ^jp the observed Young modulus is 
Yr = 2.02 X 10"^ [ksT/a^] or Yu = 2.95 x 10"^ [ksT/a^]. 



one can see an 5— shaped variation of the stress - strain relationship with initial significant elongation at vanishing 
stress due to the straightening of the membrane wrinkles (ripples) that are typical for an unperturbed membrane - cf. 
Fig. [T]d. This behavior is followed by a linear stress - strain elastic relationship where we measure the Young modulus 
Yr = 2.02 X 10~^ [ksT/a^] (or Y^ = 2.95 x 10~^ [fcsT/a^]), depending whether radial of uniaxial loading is applied. 
Eventually, for stronger stretching the elasticity of the network decreases as the anharmonicity of the bond potential 
comes into play. Moreover Fig. indicates that the destructive strain of the whole membrane is considerably less in 
the case of uniaxial stretching. 

In our work we have tried to develop the model which should serve as generic one for all kinds of 2D brittle-elastic 
networks with honeycomb orientation. We have been anxious to emphasize the common features of failure in materials 
with similar architecture but largely varying elasticity properties, e.g., from 1000 GPa graphene's Young modulus [3] 
compared to 4 x lO^'^ GPa for spectrin [T^- Putting the value of a Kuhn segment (a = 1.44 A) and taking the thermal 
energy fc^T = 4 x 10"^^ J at T = 300 K, we get from our simulation a Young modulus ~ 0.03 GPa which is ranged 
between typical values for rubber-like materials 0.01 0.1. Compared to ab initio simulations of graphene with linear 
size L = 6 which corresponds to 216 network nodes 0, our objects are about an order of magnitude larger, L = 50 
and 2500 nodes, in units of elementary cells. 
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C. MD algorithm 

As in our previous studies concerning scission kinetics of linear chains [H, [s^ and bottle-brushes [s^ we use a 
Langcvin dynamics which describes the Brownian motion of a set of interacting particles whereby the action of the 
solvent is split into slowly evolving viscous (frictional) force and a rapidly fluctuating stochastic (random) force. The 
Langcvin equation of motion is the following: 

m^(t) = (t) - m7^(t) + {t) (3) 

where m denotes the mass of the particles which is set to m = 1, iti is the velocity of particle i, i^i = M + ~^wcA)i 
is the conservative force which is a sum of all forces exerted on particle i by other particles in the system, 7 is the 
friction coefficient and Ri is the three dimensional vector of random force acting on particle i. The random force 
^i, which represents the incessant collision of the monomers with the solvent molecules, satisfies the fluctuation- 
dissipation theorem {Ria{t)Rjj3{t')) = 2'^kBT5ij5ap5{t — t') where the symbol (...) denotes an equilibrium average and 
the greek-letter subscripts refer to the x, y 01 z components. The friction coefficient 7 of the Langcvin thermostat is 
set to 7 = 0.25. Our simulation was performed in the weakly damped regime of 7 = 0.25 where effects of inertia are 
important. This value of 7 is more or less standard in Langevin MD. However, we carried out additional simulation 
in the strongly damped regime for 7 = 10. No qualitative changes were discovered except an absolute overall increase 
of the rupture times r which is natural for a more viscous environment. The integration step is 0.002 time units (t.u.) 
and the time is measured in units of rmmi/'Ti/eA/- We emphasize at this point that in our coarse-grained modeling 
the solvent is taken into account only implicitly. In this work the velocity- Verlet algorithm is used to integrate the 
equations of motion. 

Our MD simulations are carried out in the following order. First, we prepare an equilibrated membrane confor- 
mation, starting with a fully flat configuration. Fig. [1] where each bead in the network is separated by a distance 
Tmin = 1 equal to the equilibrium separation of the bond potential (Um + CAvca) [see Eq. ([1]) and ([2])]. The external 
constant force is switched on from the very beginning of the simulation. Then we start the simulation with this 
prepared conformation and let the membrane equilibrate with the applied force in the heat bath for sufficiently long 
time (w 10^ t.u.) at a temperature that is low enough so that the energy barrier for scission is high and the membrane 
stays intact. This equilibration is done in order to prepare different starting conformations for each simulation run. 
Once the equilibration is finished, the temperature is raised to the working one and wc let the membrane equilibrate 
at this temperature for roughly ~ 20 t.u. We have checked that this time interval is sufficient for equipartition and 
uniform distribution of temperature to be established throughout the membrane sheet. Then the time is set to zero 
and we continue the simulation with this well-equilibrated membrane conformation checking for scission of the bonds. 

We measure the elapsed time r until the first bond rupture occurs and repeat the above procedure for a large 
number of runs (10'^ -^ 10^), starting each time with a new equilibrated conformation so as to sample the stochastic 
nature of rupture and determine the mean (r) which we refer to as the mean first breakage time. In the course of 
simulation we also calculate properties such as the probability distribution of breaking bonds regarding their position 
in the membrane (a rupture probability histogram), the probability distribution function of the first breakage time 
W{t) (i.e., the MFBT probability distribution), the strain (extension) of the bonds with respect to the consecutive 
circle number in the membrane, as well as other quantities of interest. 

In separate runs each simulation is terminated as soon as the honeycomb sheet disintegrates into two separate parts 
whereby the time it takes to "rip-off" the sheet is termed "mean failure time (r^) and measured. In order to monitor 
the propagation of cracks, we perform also individual runs labeling breaking bonds in succession and reconstructing 
the crack trajectory which is a laborious and rather involved problem. 



D. Rupture criterion 



An important aspect of our simulation is the recombination (self-healing) of broken bonds. The constant stretching 
force acting on the monomers at the membrane edges creates a well-defined activation barrier for bond scission. Direct 
analysis of the one-bond potential with external force, J7m(£1 ~ indicates that the positions of the (metastable) 
minimum and of the barrier (or hump) r+ are given by [36| 



1 



In 



1±\/1-/' 



(4) 



5 



where the dimensionless force / = 2//aej\/. For the range of tensile forces used in the present study one has typically 
r_|_ « Srniin. The activation energy (barrier height) for single bond scission is itself given by [36j 
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One can easily verify that E't, decreases with /. Since a bond may get stretched beyond the energy barrier and 
nonetheless shrink back again, i.e. recombine, in our numeric experiments wc use a sufficiently large value for 
critical extension of the bonds, r/j = 5r,„i,i, which is defined as a threshold to a broken state of the bond. This 
conventions is based on our checks that the probability for recombination (self-healing) of bonds, stretched beyond 
r/i, is vanishingly small, as demonstrated below. In our model wc deal with Eb/{kBT) ~ 20 which at 300 K and 
bond length r,„i„ = 0.14 nm corresponds to ultimate tensile stress ^ 0.6 GPa. This is a reasonable value for our 
membrane which is considerably softer than grapliene with ~ 100 GPa [3| and is ranged between typical values for 
rubber materials 0.03 ^ 14 GPa. 



III. MD-RESULTS 



We examine the scission of bonds between neighboring nodes in the network sheet with honeycomb topology, 
assuming thermal activation as a driving mechanism in agreement with early experimental work by Brenner |37| and 
Zhurkov [H]. In Fig. [3] we show a series of representative snapshots of a membrane of size i = 10 with N = 600 
monomers taken at different time moments during the process of decomposition. Typically, the first bonds that break 
are observed to belong to the last (even) most remote circle as, for example, aX t k, ITlt.u. in Fig. |3l As mentioned 
above, these are the radially oriented bonds which belong to concentric circles of even number. Gradually a line of 
edge beads is then severed from the rest of the membrane and a crack is formed which propagates into the bulk until 
eventually a piece of the network sheet is ripped off, as in Fig. [3] at t w 370i.u. As we shall see below, this mechanism 
of membrane failure, whereby an initial crack is formed parallel to the edge monomers, yet perpendicular to the tensile 
force, dominates largely the process of disintegration under constant tensile force. The process is, therefore, mainly 
described by two characteristic times, r and t^, which mark the occurrence of the first scission of a bond (MFBT) 
and that of the eventual breakdown of the flake into two distinct parts. 



A. Bond recombination 



As mentioned in Section Hill throughout in our studies of the brittle sheet breakdown wc use a threshold for critical 
bond stretching (rupture criterion) ~ 5r,„i„. In the right inset of Fig. 2] we display the function Qh{h), which 
represents the probability distribution of bond stretching h beyond the hump position r+ , given that a subsequent 
recombination has taken place. To this end one monitors for 10** integration steps the length of each bond once the 
bond expands beyond r+ and stores its maximal expansion, h, provided such a bond contracts again to r < r+. 
Then Qh{h) is computed as the fraction of extensions to h over the total number of recombination events. For each 
bond recombination one measures also the distribution of the respective self-healing times, /\(t), which is shown 
in Fig. 2] too. Both distributions are characterized by exponentially fast decaying tails, indicating that successful 
recombinations are possible after very short time interval « l.Zt.u., and the possible stretching of a bond in such 
cases is minimal - about 0.19-^0.5 beyond the energy barrier position at r+ w 2.96, that is, significantly shorter than 
r/j K, 5. We also find that recombination of bonds takes place seldom (roughly 1.5% over 5 • lO"' runs of average length 
K, Ayit.u. for a membrane composed of = 600 beads). Yet as indicated below, allowing for self-healing events may 
significantly change the observed kinetics of membrane destruction. The left inset in Fig. 2] indicates that self-healing 
of bonds happens most frequently at the membrane periphery, C = 19, where bond stretching occurs most frequently. 



B. Mean First Brealcage Time 



These conclusions, based on visual evidence from snapshots taken in the course of membrane decomposition, are 
corroborated in Fig. [5ji where we show the probability distribution of a first rupture for all bonds in the honeycomb 
membrane flake as a 3D plot. It is seen that the scission rate is localized in the outer- most circle of radial bonds 
whereas bonds in the inner part of the membrane practically hardly break. Note that this is not a trivial effect since 
tension is distributed uniformly over all bonds in the equilibrated membrane so there is no additional propagation of 
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FIG. 3: Snapshots illustrate the process of bond breakage (crack generation) in different time moments for a membrane with 
= 600 particles subject to external tensional force / = 0.15 at T = 0.05 and 7 = 0.25. The force is applied to periphery 
monomers only and stretches the network perpendicular to its original edges. 



the tension front from the rim towards the center. Fig. \5jp also indicates a qualitative change in the rupture PDF 
when self-healing is not allowed for (by reducing the threshold position to that of the energy barrier - = 3.1) in 
contrast to results where self-healing was fully taken into account - rh = 5. Moreover, a closer inspection the new 
Fig. [SJd indicates that scission of bonds with no self-healing takes place almost uniformly throughout the membrane 
while with self-healing it is concentrated only at the membrane periphery. 

One can try to relate this finding to the distribution of strain within the network as shown in Fig. and sampled 
for several strengths of the external stretching force /. In the case of strongest pulling, / = 0.15, the variation of 
the mean-squared bond length {P) with distance from the membrane center (i.e., with consecutive circle number C) 
displays a well expressed saw-tooth behavior whereby the peaks correspond to bonds with radial rather than tangential 
orientation (odd C). The alternation of strongly / weakly stretched bonds modulates the overall gradual increase of 
the mean bond length with growing distance from the center. Evidently, the amplitude of the mean-squared bond 
length attains a pronounced maximum on the last circle of radially oriented network bonds. This distribution of strain 
is found to persist down to vanishing tensile force / = - Fig. [5^. The distribution of first scission events is clearly 
seen in Fig. [6)3 where we show it for several strengths of /. Evidently, with growing value of / bonds happen to break 
also deeper inside the membrane although such events remain much less probable. 
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FIG. 4: Probability distribution Ph{t) of maximal times (full circles), and Qh{h) of maximal bond lengths h (circles, right 
panel of inset) before a recombination event in the stretched membrane with A'" = 600, T = 0.1, 7 = 0.25 takes places. The 
exponential tail of Ph{t) is fitted by blue line. The exponential decay of Qh(h) is given by red line. The left panel of the inset 
shows the healing probability Rh vs. circle number C. The healing events under applied stress occur roughly 10 times less 
frequently than for / = 0. 
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FIG. 5: (a) Rupture probability histogram of flexible hexagonal membrane subjected to external tensile stress / = 0.15. (b) 
Scission probability histogram vs consecutive circle number for membrane pulled with / = 0.125 displayed for different rupture 
thresholds rh as indicated. Real rupture events [rn ~ 5.0) are concentrated at the periphery whereas fictitious ones (rh = 3.1) 
are distributed more uniformly all over the membrane. Here — 600, T — 0.05 and 7 = 0.25. 



The variation of the MFBT r with system size N (i.e., with the number of monomers in the membrane N ~ 6L^ 
where L denotes the hnear size of the flake) is shown in Fig. [71 For sufficiently large membranes one observes a power 
law decline of the MFBT, t oc with an exponent /3 « 0.5 ± 0.03 for the tensile forces studied.. If thermally 

activated bonds break independently from one another and entirely at random, then the MFBT t measures the 
interval before any of the available intact bonds undergoes scission, that is, either the first bond breaks, or the second 
one, and so on which, at constant rate of scission, would reduce the MFBT r oc 1/N as observed for instance in the 
case of thermal degradation of a linear polymer chains (3^ . A more comprehensively this simple result can be derived 
by means of the classical theory of WeibuU. In the present system of a honeycomb membrane the bonds that undergo 
rupture are nearly all located at the rim of the flake and their number is proportional to L so that with /3 « 0.5, 
cf. Fig. [7] and N oc L'^, one obtains eventually the important result t oc 1/L. This observation is in agreement with 
recent results of Grant et al. [2^ who studied the nucleation of cracks in a brittle 2D-sheet. We should like to point 



8 



V. 






0.15 


f = 


0.125 


o-of — 






^ \ , I , \ , \ , \ ^ 

10 20 30 40 50 

Consecutive circle number 




FIG. 6: (a) Variation of the mean squared bond length {l^} with consecutive circle number in a membrane with A'^ = 5400 
beads subjected to different strengths of the external force / (as indicated in the legend) . (b) Probability distribution of the 
first bond scission event vs circle number in a membrane with A'^ = 600 beads at different strengths of the external force / as 
indicated. For a force / < 0.15 the bonds from the last two outer circles (#18 and #19) in the membrane have the highest 
rupture probability. With increased strength of the pulling force / > 0.175 the bonds which are located in the circles #18 and 
#16 attain the highest rupture probability. Parameters of a heat bath are T = 0.05 and 7 = 0.25. 



out at this place that without self-healing, c.f., Fig. [Sj), rupture time goes as (r) oc N~^, with /3 w 1 (not plotted 
here) in contrast to the observed exponent /3 « 0.5. 

One should mention here an interesting analogy between the observed power-law dependence of the MFBT time on 
system size and the power-law decrease of life-time with system size in thermally activated breakdown of fiber bundles 
(30I [sij . While both in our honey-comb network as well as in the Fiber Bundle Model (FBM) the failure mechanism 
is related to redistribution of load on neighboring bonds (fibres) upon single rupture, bonds in our membrane are 
subject to a single scission threshold whereas in the FBM there is a random distribution of tensile strengths. As a 
result, one finds a single value of /3 w 0.5 and an Arhenian dependence of r on temperature T in our elastic-brittle 
honey-comb network (see below) while the exponent /3 depends on the external load / and on T giving rise to a 
non- Arhenian t vs. T relationship. 

Note that the decline of MFBT r in a topologically connected brittle system is by no means a trivial one. In 
a recent study [1^, using Molecular Dynamics (MD) simulation of a single anharmonic polymer chain subject to 
constant external tensile force, we found a rather complex interplay between the polymer chain dynamics and the 
resulting bond rupture probability distribution along the chain backbone. In a breakable chain (rather than 2D 
network) it was observed that the corresponding power /3 — ;> as iV — ^ 00. A major role in this was attributed 
to nonlinear excitations as the possible origin for the observed increasing insensitivity of rupture time with respect 
to polymer length as the pulling force grows. One may thus conclude that nonlinear effects in bond scission are 
suppressed in 2D honeycomb networks. 

One can also see from the inset in Fig. [TJi that the MFBT r decreases rapidly with growing stress /, that is, the 
energy barrier for rupture declines with / in agreement with Eq. ([5]) and Zhurkov's experiments [38|. The probability 
distribution of MFBT W{t) is shown in Fig. [7l3. It is well described by a Poisson probability distribution function 
Wit) = 5.57 • lO^H^/^ exp(-t/291.85). 



C. Cracks and Mean Failure Time 



The variation of r^, the mean failure time of the membrane with system size N, shown in Fig. [5Ji, displays also 
a power-law dependence on system size N, (t,.) cx N~'^, whereby (p undergoes a cross-over to a lower value beyond 
roughly N > 300. However, has different physical meaning. Following Pomeau [13, the failure time can be 
approximately identified with the nucleation of a crack of critical size Ic given by Griffith's critical condition [l^, [ioj 
assuming that crack propagation is much faster than the nucleation time. For a 2£)-geometry consisting of a fiat 
brittle sheet with a crack perpendicular to the direction of stress, the potential energy per unit thickness of the sheet 
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FIG. 7: (a) Mean first breakage time (r) vs. number of particles A*' in tiie membrane pulled with different tensile stress / 
as indicated. Symbols represent simulation data whereas solid lines stand for fitting functions (r) ~ . The inset shows 
force-dependent (r) for a membrane composed of A'^ = 294 beads, (b) MFBT probability distribution W{t) for the first scission 
of a bond in a flake with A'' = 600 particles at stress / = 0.15. Symbols denote result of simulation and full line stands for the 
fitting function W{t) oc t" exp(-t/r) with n = 1/3 and r = 291 -8.5/;. 7/,. Parameters of heat bath are T = 0.0.5 and ^ = 0.25. 
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FIG. 8: (a) Mean failure time (jr) (time needed to split membrane into two pieces) vs. number of particles in the membrane for 
two values of the external pulling force / at T = 0.05 and 7 = 0.25. Symbols denote simulation results and solid line represents 
power law fitting function {jr) ~ N"'^. The inset shows PDF of number of particles in the moment of splitting for a membrane 
composed of A'' = 150 beads, (b) Failure time (rr) vs / in the case of A'^ = 294. 



reads U ~ — + 2eZ + Uq where Y is the Young modulus, e is the surface energy needed to form a crack of length 
Z, and Uq is the elastic energy in the absence of stress (/ = 0). This energy reaches a maximum for a critical crack 
length Ic = ^ji beyond which no stable state exists except the separation of the sheet into two broken pieces. Thus, 

with a crack nucleation barrier AC/ — (in 3D AC/ oc Z^^'), the failure (rip-off) time = tq exp(AC/o/fcBr) as 

found in experiments with bidimensional micro crystals by Pauchard and Meunier (4]| and in gels by Bonn et al. (4^ . 
In Fig. [8]d we present the variation of for membrane failure with stress / in good agreement with the expected 
relationship AC/ cx /~^. In addition, we show the variation of Tr with temperature, see inset in Fig.|8l), which is found 
to follow a well expressed Arhcnian relationship with inverse temperature, in agreement with earlier studies [29l |40|. 

The end of the sheet rupturing process is marked as a rule by disintegration into two pieces of different size so it is 
interesting to asses the size distribution of such fragments upon failure. In the inset in Fig. |8K we show a probability 
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distribution S{n) of the sizes of of both fragments upon membrane rip-off. In a membrane composed of N beads one 
observes a sharp bimodal distribution with narrow peaks at sizes A^i « 10 and N2 « 140. Evidently, for the adopted 
nearly radial direction - cf. Fig. [2^ - of the applied tensile force one always finds a pair of one small and another very 
large fragment. 
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FIG. 9: Overview of observed cracks in a honeycomb membrane composed of A'' = 600 particles for different orientation of 
the applied external pulling force. Green arrows show indicate the orientation of the applied force (/ = 0.15): (a), (b), (d) - 
uniaxial, (c) - biaxial, (e), (f) - slanted. Parameters of a heat bath are T — 0.05 and 7 = 0.25. The typical cracks are marked 
in color on the geometrically undistorted arrangement of network nodes for better visibility. 



One can readily verify from the typical topology of the observed cracks in the membrane, presented in Fig. |9l that 




FIG. 10: (a) Crack propagation velocity (number of broken bonds per unit time) for a membrane with N — 600 beads at 
different strength of the external force / as indicated, (b) Three different realizations of cracks at applied force / = 0.14. The 
inset shows a variation of the mean crack propagation velocity with /. Here T — 0.05 and 7 = 0.25. 



(i) cracks emerge as a rule perpendicular to the direction of applied stress, and (ii) it is almost always the first row 
of nodes to which the tensile force is immediately applied that gets ripped off upon failure. Cracks that break the 
network sheet in the middle occur very seldom, in compliance with the sampled distribution of fragment sizes, S(n) 
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in the inset of Fig. |8^. One would, therefore, predict a breakup of a protective cover spanned on the orifice of tube 
like the one shown in Fig. [5] to proceed immediately at the fixed orbicular boundary where the tensile force applies 
to the network. It is interesting to note that the geometry of cracks in the membrane shown in Fig. |9] appears very 
similar to the one observed in drying induced cracking of thin layers of materials subject to structural disorder (32| . 

The emerging cracks are expected to propagate with speed that increases as the strength of the external force is 
increased as the inset in Fig. llOb , indicates. In fact, in Fig. llOb one observes typical curves comprising a series of short 
intervals with steep growth of the number of broken bonds per unit time and longer horizontal 'terraces' preceding 
the nucleation of a new crack. Even though the data, presented in Fig. llOb . is not averaged over many realizations, 
and, as Fig. [TOb suggests, individual realizations of propagating cracks may strongly differ even at the same stress 
/, a general increase of the propagation velocity with growing external force / - see inset - can be unambiguously 
detected, in agreement with earlier observations [20| . 

For our model membrane with computed Young modulus Y w 0.02 we get for the Rayleigh wave speed cr ~ 0.14. 
Thus for most of the applied tensile stress values we observe crack propagation at speed less than cr - inset in Fig.lTOb. 
As argued by [i^ propagation speed cannot exceed cr because crack splits off into multiple cracks before reaching cr. 
In contrast, Abraham and Gao in Ref. [31 have reported on cracks that can travel faster than the Rayleigh speed. 
Thus, our rough estimates (inset in Fig. llOp agree well with data from literature. Converting our results to proper 
metric units, with bond length cr ss 0.144 nm and energy « 2QkBT which yields 1 MD t.u. w 10^^^ s, we estimate 
the typical crack propagation speed Vc ~ 50 m/s. Note that mean crack speed for natural latex rubber was given as 
56 m/s m. 

IV. SUMMARY 

In the present work we have studied the bond rupture and ensuing fracture of a honeycomb brittle membrane subject 
to uniform radially applied external stretching forces for different values of force /, temperature T, and membrane size 
N. The most important conclusions that can be drawn from our Molecular Dynamics simulation can be summarized 
as follows: 

• bonds scission in hexagonal 2D sheets with honeycomb structure of the underlying network under subjected to 
external pulling perpendicular to flake's edges take place overwhelmingly at the sheet periphery 

• The Mean First Breakage Time of breaking bonds depends on membrane size A'^ as a power law, t oc with 
P « 0.50 ±0.03. 

• The failure time until a brittle sheet disintegrates into pieces follows a power law too, r a A^~"^W), and an 
exponential decay oc exp(const//^) upon increasing strength of the pulling force, in agreement with Griffith's 
criterion for failure. 

• cracks emerge in the vicinity of membrane edges and typically propagate parallel to the edges, splitting the 
sheet in two pieces of size ratio of « 7%. 

• crack propagation speed is observed to increase rapidly with tensile force 

We believe that these findings can be seen as generic also for 2D network brittle sheets of different geometry 
(hexagonal lattices, or quadratic lattices with second nearest-neighbor bonding) where similar interplay between 
elastic and fracture behavior is expected to take place. It is clear, however, that more investigations are needed before 
a full understanding of fracture in such systems is achieved. 
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